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Abstract 

A numerical hydrodynamic study of femtosecond laser ablation is presented. A detailed analysis 

o. 

O , of material decomposition is performed using a thermodynamically complete equation of state with 

1*3 i 

separate stable and metastable phase states and phase boundaries. The lifetime of the metastable 
liquid state is estimated based on the classical theory of homogeneous nucleation. In addition, 

43 ; 

mechanical fragmentation of the target material is controlled based on available criteria. As a result, 
several ablation mechanisms are observed. A major fraction of the ablated material, however, is 
found to originate from the metastable liquid region, which is decomposed either thermally in 
the vicinity of the critical point into a liquid-gas mixture, or mechanically at high strain rate 
. and negative pressure into liquid droplets and chunks. The calculation results explain available 

[ — . ' experimental findings. 

o ' 
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I. INTRODUCTION 



Femtosecond laser interactions have attracted a particular interest of researchers during 
the last decade. It is not surprising since these new ultra-short laser pulses provide a unique 
opportunity for the fast development of numerous very promising and exciting applications, 
such as laser machining, micro and nano-structuring, laser-induced plasma spectroscopy, 
nanoparticle synthesis in vacuum, in gas or in a liquid solution, medical imaging, laser 
surgery, etc. l], [2], [J Many experimental studies were performed demonstrating the unique 
possibilities of femtosecond laser systems. At the same time, several theoretical 



investigations were focused on the better understanding o 



ultra-short laser interactions. 
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2l| The 



classical approach typically used for these interactions is based on a two-temperature model. 
[13I In addition, several hydrodynamic simulations were carried out to describe the target 
material motion, [h], [ll, 14, Q] The formation of a thermal wave, pressure gain and 
shock wave propagation have been observed in these calculations. Molecular dynamics (MD) 
simulations were furthermore performed that provided insights into the mechanisms, such as 
phase explosion, fragmentation, evaporation, and mechanical spallation. [l7, Q, 
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In spite of a significant effort aimed at the understanding of the ablation mechanisms, sev- 
eral issues stay unclear. For example, theoretically predicted ablation depth was frequently 
underestimated with respect to the typical experimental values. This result is an indication 
that not all the ablation mechanisms were taken into account. In addition, numerous terms 
used to describe these mechanisms are confusing or misused. For instance, when explosive 
boiling term is used, which reflects the overheated (metastable) liquid decomposition into 
a liquid-gas mixture, it is often associated with material heating. The metastable liquid 
state, however, can be also achieved during fast material expansion and cooling. The proper 
term for this process is unclear. In this case, critical-point phase separation was proposed 
in Ref. UM, but only to describe material decomposition in the vicinity of the critical point. 
In addition, up to now, it was difficult to handle this process correctly in the framework 
of hydrodynamic models, because most of the available equations of state (EOS) did not 
describe the metastable regions. Moreover no kinetics of gas bubble formation or/and con- 
densation was previously introduced to determine the lifetime of these states. Besides, there 
is still a lack of information about the photomechanical failure of the material at extremely 
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high strain rates (up to 10 9 s 1 ), which can be achieved. Nevertheless, even these high 



231 ] because of their short 



strain rates do not always satisfy the fragmentation criteria [2 
duration (tens of picosecond). All these issues need more analysis. 

II. MODEL 

To bring more light on the mechanisms of material decomposition in femtosecond laser ab- 
lation of metals, we performed a detailed hydrodynamic modeling. Remind, that compared 
to atomistic techniques, hydrodynamic models allow calculations for much larger systems 
and take much shorter computer time. Our numerical model is based on the solution of a sys- 
tem of Eulerian hydrodynamic equations by a high-order multi-material Godunov's method. 



24| The equations were extended to the case of one-fluid two-temperature hydrodynamics 



with laser energy absorption source, electron heat conductivity and energy exchange be- 



tween electrons and heavy particles (atoms, ions, nuclei). 25[ A Gaussian temporal profile 



is used to simulate the laser energy deposition. The heat conductivity of electrons is calcu- 



lated according to model. 26] Energy exchange between electrons and heavy particles, the 



reflectivity coefficient R and the optical penetration depth \ opt are derived by means of the 
wide-range frequency of electron-phonon/ion collisions, [l^] 

For completeness of our model, we use a semiempirical multi-phase EOS for aluminum 
with separate description of subsystems of heavy particles and electrons. The specific 
Helmholtz free energy has a form J-(p, Tj,T e ) = J~i(p, Tj) + jF e (p, T e ), composed of two 
parts, which describe the contributions of heavy particles and electrons, respectively. Here 
p is the material density, Tj and T e are temperatures of heavy particles and electrons. The 
first item !Fi(p, Tj) = T c {p) + T a (p, Tj), in turn, consists of the electron-ion interaction term 
T c (calculated at Tj = T e = K) and the contribution of thermal motion of heavy particles 
T a - Analytical form of Ti has different expressions for the solid as well as for both 
liquid and gas phases J-f . |27] Using these thermodynamic functions, the solid, liquid, and 
gas phases equlibrium boundaries are determined from the equality conditions for the tem- 
perature Tj, pressure Pj = p 2 (dJ r i /dp)T i , and Gibbs potential Gj = Ti + Pi/ p of each phases 



pairs. 



The tables of thermodynamic parameters are calculated ta king into consideration 



the information about phase transitions and metastable regions. 29 
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electrons in metal T e has a finite-temperature ideal Fermi-gas form. 28] 



The free energy of 
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The model of our EOS contains about 40 adjustable constants so that the EOS meets the 
following requirements: (i) to describe experimental results on compression and expansion 
for a wide range of densities and temperatures including data on critical and triple points; (ii) 
to contain separate information about electron and ion/lattice subsystems; (iii) to represent 
changes of thermodynamic parameters during phase transitions. To account for kinetic 
processes, an estimation of the realistic lifetime of metastable liquid state is introduced. 
Using this lifetime-based information we switch between two different modifications of the 
EOS: (i) with metastable states; and (ii) without metastable states. 

In our model, when the liquid branch of the binodal curve was crossed and the matter 
turned into metastable state, we include a particular treatment for each of the following 
two competitive effects: (i) for the thermal decomposition, a criterion of the metastable 
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33| (ii) for 



liquid lifetime is used based on the theory of homogeneous nucleation- 
the mechanical fragmentation, a failure criterion of Grady is applied. [22j 

In the first case, we estimate the metastable liquid lifetime as r nuc i = 
(CnV)~ 1 exp(W/kBTi), where C = 10 10 s _1 is the kinetic coefficient, [33J] n is the con- 
centration, V is the volume under consideration (the volume of numerical cell in our case), 
W = 16/T(j 3 /3AP 2 is the work required on a gas bubble formation in a liquid phase, ks is 
the Boltzmann constant, AP is the difference between the saturated vapor pressure, known 
from EOS, and the pressure of substance. The temperature dependence of the surface ten- 
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34. Here T r = 6595 K is 



sion a is described in the form a = a (l — Tj/T r ) , see Ref. 
the temperature in critical point (CP) for Al, |27j ctq = 860 g/s 2 is the surface tension of Al 
at melting temperature. 34J In this case, as soon as the lifetime r nuc \ in the volume V is 



expired, the metastable one-phase state separates into two-phase stable liquid-gas mixture. 
This phase separation process is accompanied by abrupt changes of thermodynamic param- 
eters such as pressure, temperature, sound speed, compressibility, heat capacity, etc. The 
EOS with metastable phase states is therefore no more relevant in this volume, so that we 
continue to calculate the thermodynamic properties by using the EOS without metastable 
states. 

I — I 

To account for the second effect, a fragmentation criterion [22j is used for the liquid phase 
with the spall strength P spa i = (6p 2 c 3 cre) 1//3 and the time to fracture r spa i = (6a/p£: 2 c 3 ) 1 / 3 , 
where e is the strain rate and c is the sound speed. When the pressure of matter drops 
below the negative value —P spa i and duration of this event is longer than the time T spa i, 
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a criterion of fragmentation is satisfied and voids and new free surfaces can appear in the 
substance. To describe this phenomenon numerically we let a "broken" substance to shrink 
back until the pressure comes to zero value. The difference between old and new volumes we 
compensate by introducing void fraction into a numerical cell. Both thermal and mechanical 
criteria described above are used simultaneously and each of them can prevail in a given 
computational cell depending on the substance location on the phase diagram. 



III. RESULTS 



For the analysis of the target material evolution after the laser irradiation, it is convenient 
to use a phase diagram in the temperature-density plane (T,p) given by the EOS (Fig. [[]). 
The binodal curve corresponds to the liquid-vapor coexistence line (saturated vapor curve). 
The spinodal curve (dashed) gives a limit of thermodynamic stability of the liquid and gas 
phases. The regions between the binodal and the spinodal curves show the metastable states 
of matter: the superheated liquid (1) and the supercooled vapor (g) states. 

Initially, the aluminum target is in a solid state with p = 2.71 g/cm 3 , Tj = T e = 300 K, 
and is subjected to a laser pulse with ffuence F = 5 J/cm 2 , pulse width tl = 100 fs and 
wavelength Ar = 800 nm. The reflectivity coefficient and the optical penetration depth are 
calculeted [14J and their values are R = 0.92 and X opt = 13 nm, respectively. The conduction- 
band electrons absorb laser energy in a skin-layer (a very thin region of the order of the 
optical penetration depth A op t). Then, these hot electrons transmit absorbed energy into 
the bulk owing to heat conduction and warm up the lattice due to electron-phonon collisions. 
The time of electron-lattice relaxation (equilibration of temperatures of electron and heavy 
particle subsystems) is of the order of several picoseconds for aluminum and thus the lattice 
pressure has time to grow up before the heat energy of electrons dissipates into the bulk of 
the target. Formation of a narrow zone with a sharp pressure profile results in generation of 
intensive shock wave, which moves into the bulk, and rarefaction wave that causes a spread 
of the matter out of the target. The process of the shock wave propagation is clearly visible 
as the consecutive right-hand deviations of the 30, 50, 80, 130 nm-trajectories in Fig. [IJ 
Process of compression of different target layers changes into expansion one (decreasing of 
the density along trajectories, Fig. CD). 

The maximal temperature of heavy particles in the vicinity of the target surface (depth 
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FIG. 1: Phase diagram of aluminum and time evolution of density and temperature of heavy 
particles for different layers of the target. Dashed curve: the spinodal; g: stable gas; 1: stable 
liquid; s: stable solid; 1+s: stable melting; g+1: liquid-gas mixture; g+s: solid-gas mixture; (g): 
metastable gas; (1): met ast able liquid; (1+s): metastable melting; (s): metastable solid. The phase 
trajectories 1-7 correspond to depths of 5, 15, 20, 30, 50, 80, 130 nm from the initial target surface, 
respectively. Arrows along the trajectories show the flow of time. Here, laser pulse parameters are 
t l = 100 fs, X L = 800 nm, and F = 5 J/cm 2 . 

< 15 nm) is high enough (up to Tj ~ 30 kK), so that the phase trajectories from this 
layer go above the CP and the target material is directly transformed into the gas phase 
(atomization, 5 and 15 nm trajectories in Fig. [p. Then, these trajectories cross the gas 
branch of the binodal and penetrate into the supercooled vapor region, where condensation 
starts and liquid-gas mixture forms. Note, however, that this layer represents a very small 
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fraction of the ablated material, so that the condensation degree is insignificant. 

Next layer in the target (depth from 20 to 30 nm) is first transformed into a metastable 
liquid state. This layer is heated to high temperatures (up to Tj ~ 25 kK), so that the 
corresponding phase trajectories cross the liquid branch of binodal in the vicinity of the 
CP and enter the metastable liquid region. Under these conditions, the target material is 
thermodynamically unstable. The lifetime of this state is estimated as described above. 
Thermodynamic instabilities are known to occur near the CP (particularly, at 0.9T C < Tj < 
T c ) leading to a rapid decomposition (several picoseconds) of the matter into a liquid- 



gas mixture. This process is similar to the phase explosion, 



17 



18] though occurs during 



simultaneous material expansion and cooling. The other term used for this mechanism is 
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critical- point phase separation proposed in Ref. 

This thermal mechanism, however, concerns only quite a thin slice (~ 10 nm thickness), 
which trajectories enter the metastable liquid region near the critical point. In fact, the 
melted layer is much thicker (~ 200 nm at F = 5 J/cm 2 ). The rest of the trajectories (5-7 
in Fig. CP, which originate from the deeper lying melted layers of the target, also enter the 
metastable liquid region, but at temperatures much smaller than the critical temperature T c . 
In this case, the matter can stay in the metastable state much longer. As a result, mechanical 
effects start to dominate over the thermal ones and material fragmentation (cavitation in 
liquid) occurs under the action of a tensile pressure wave. In other words, the time to 
mechanical fracture is shorter than the nucleation time {jnuci < T spa i). As soon as the time 
to fracture is exceeded, the material is decomposed into droplets and chunks (separated by 
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concerns the 



vacuum). This mechanical decomposition, described also in Ref. 1?], 
major part of the ablated material in our simulations. 

The dynamics of the target surface layers decomposition is presented in Fig. [2j Initially, 
the solid-vacuum interface is located at x = nm and target occupies x > nm half 
space. During the collisions of heated electrons with the lattice its temperature and pressure 
grows, and thus shock wave forms and propagates into the target (pressure peaks in Fig. [2j) . 
Pressure release on the free surface of the target results in material expansion in the opposite 
direction through a rarefaction wave. At a delay t ~ 10 ps after the beginning of the laser 
pulse, the density profile is still continuous [Fig. [2]^a)]. At t ~ 20 ps, the trajectories 
originated from the skin-layer reach the binodal curve and penetrate into the metastable 
liquid region, where thermal decomposition of substance starts [Fig. [2(b)] . This process is 
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FIG. 2: Contour plots of density (solid) and pressure of heavy particles (dashed) for different time 
delays after irradiation: (a) — 10 ps, (b) — 20 ps, (c) — 30 ps and (d) — 80 ps. Here, 1 — is thermal 
decomposition zone, 2 — is mechanical fragmentation zone. Laser parameters are the same as in 

Fig.m 

completed by t ~ 30 ps [Fig. [2(c)] . Then, only mechanical effects, caused by the tensile 
stress (negative pressure), are involved into the process of decomposition [Fig. EJ^d)]. By 
the time t ~ 80 ps the decomposition of the liquid phase is completed, because the tensile 
pressure wave intensity drops owing to the energy expenditure on fragmentation (compare 
minimal pressures in Fig. [2](c)-(d) at x ~ 40 and x ~ 260 nm, respectively). Completion of 
the mechanical decomposition corresponds to the final vibrations of the trajectories 5-7 in 
Fig.[T]around isobar with zero pressure, similar to recent MD analysis [21] and hydrodynamic 
simulations. [g] 

The fractions of the target material ablated due to the described mechanisms depend 
both on the material and on laser parameters. At very small fluences (F < 0.25 J/cm 2 for 
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aluminum) only melting occurs. When laser fluence is slightly above this value, material 



fragmentation takes place, in agreement with the previous MD simulations. [12j, [20J If we 
further increase laser fluence, both mechanical and thermal mechanisms play a role. At 
larger F, all three mechanisms (i) atomization; (ii) thermal decomposition (critical-point 
phase separation); and (iii) mechanical decomposition (dynamic fragmentation) occur in 
different layers of the target. The third mechanism, however, is found to be dominant in our 
simulations, accounting about 80% of ablated mass at all considered laser fluences (from 0.1 
to 5 J/cm 2 ). 

Finally, we propose several predictions from our simulations and compare our results 
with the available experimental findings 5|, |7| and simulations. [loT ] We use a single- shot 
experimental data for the numerical model verification, though oxidation of aluminum target 
surface and usage of different methods of measurement result in essential discrepancy in 
single-shot laboratory findings, see Fig. [3j It turns out, that the consideration of only critical- 

gives rise to a significant underestimation 
of the ablation depth (several times) with respect to the experimental values (Fig. [3]). This 
underestimation can be explained by the fact that only 10-20% of the target material is 
ablated due to these critical effects, whereas most of the ablated material is ejected due to 
the mechanical fragmentation of the liquid phase. These results confirm that the decomposed 
melted zone contributes strongly into the estimation of the ablated depth. 



IV. CONCLUSION 



In conclusion, we have presented the results of the hydrodynamic calculations of femtosec- 
ond laser ablation with a thermodynamically complete equation of state and with interplay 
between both kinetic lifetime of metastable liquid state and mechanical time to fracture. 
Ablation mechanisms have been investigated as a function of laser fluence. A correlation 
between the melted and ablated depths has been revealed. Mechanical decomposition of 
a metastable liquid layer has been shown to be the dominant ablation mechanism. In the 
present work we have considered only aluminum ablation in one dimension. Similar anal- 
ysis can be performed for different materials. In addition, an extension of the model to a 
two-dimensional case with multi-layer materials is underway. 
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FIG. 3: Present simulation: melted (A, solid line) and ablated (a) depths as a function of laser 
fluence for aluminum. Experimental data are taken from Ref. |5j — (■) and Ref. |7| — (D)- The 
errors of the experimental data were reported to be of the order of 5-10%. Results of simulation 
(*) and Ref. ▼ ) are also shown. 
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